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ABSTRACT 

We conduct a theoretical study on the ejection of runaway massive stars from R136 
— the central massive, star-burst cluster in the 30 Doradus complex of the Large Mag- 
ellanic Cloud. Specifically, we investigate the possibility of the very massive star (VMS) 
VFTS 682 being a runaway member of R136. Recent observations of the above VMS, by 
virtue of its isolated location and its moderate peculiar motion, have raised the funda- 
mental question whether isolated massive star formation is indeed possible. We perform 
the first realistic N-body computations of fully mass-segregated R136-type star clusters 
in which all the massive stars are in primordial binary systems. These calculations 
confirm that the dynamical ejection of a VMS from a R136-like cluster, with kinematic 
properties similar to those of VFTS 682, is common. Hence the conjecture of isolated 
massive star formation is unnecessary to account for this VMS. Our results are also 
quite consistent with the ejection of 30 Dor 016, another suspected runaway VMS from 
R136. We further note that during the clusters' evolution, mergers of massive binaries 
produce a few single stars per cluster with masses significantly exceeding the canonical 
upper-limit of ISOM©. The observations of such single super-canonical stars in R136, 
therefore, do not imply an IMF with an upper limit greatly exceeding the accepted 
canonical ISOM© limit, as has been suggested recently, and they are consistent with the 
canonical upper limit. 

Subject headings: stellar dynamics — methods: N-body simulations — stars: individual 
(VFTS 682, 30 Dor 016) — galaxies: star clusters — open clusters and associations: 
individual (R136) 



1. Introduction 

Fast-moving massive field stars in our Galaxy and in external galaxies, that apparently are 
not associated with any nearby stellar assemblies, have been a focal topic for the past two decades. 
Such runaway OB stars are widely thought to be members of star clusters that are ejected following 
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dynamical encounters in their parent clusters. If the proper motions of these massive stars can be 
reliably estimated, the v ast majority of thern can th en be traced back to their parent clusters. Using 
this kinematic method, ISchilbach &: Rosed (|2008l ) have demonstrated that most of the Galactic 
OB runaway stars indeed originated from star clusters. Another way to determine the parents of 
runaway massive st a.rs is to image their bow shocks , if present, which has been possible with modern 

A detectable bow shock is generated if a star 
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infrared telescopes (jGvaramadze et al.l 12010 . 
moves supersonically (typically > 10 km s~^) through a dense enough medium. The geometry of a 
bow shock allows one to d etermine the direction of motion of the shock-generating star and hence 
its possible parent cluster (jGvaramadze et al.ll20ld ). Bow shocks are particularly useful when the 
proper motions of the runaway OB stars are unavailable or are poorly known which is often the case 
for distant stars. However, even if every OB field star has been ejected from a young cluster, 1-4% of 
the OB stars cannot be traced back t o their parent clusters because of the two-step ejection process 
(jPfiamm-Altenburg &: Kroupall2010l ): a massive binary is first ejected from its cluster and when 
the primary explodes as a supernova, the secondary OB star continues on a diverted trajectory. 

High- velocity single and binary stars are launched fr om a star clu ster due to super-elastic single- 
star — binary and binary-binary dynamical encounters (|Heggidll975l ) that occur predominantly in 
the cluster's dense core. Therefore, the s pectrum of the ejected bodies is governed by the pr operties 



of the primordial binaries in the cluster ([Leonard Duncan 



1988 



Clarke fc Pringlelll992l ) and its 



stellar initial mass function (IMF) and hence serves as a useful tracer of these fundamentally 
important cluster properties. As the stars and the binaries are mostly ejected fro m the cluster 's 
core, where the most massive members remain segregated (mass stratification; see ISpitzed 119871 ) . 
the runaway population contains a large fraction of OB single stars/binaries. The ejection of the 
massive members also alters the stellar mass function (MF) of the bound cluster as it evolves. 

A particular runaway very massive star (VMS) that has recently drawn significant attention 
is VLT-FLAMES Tarantula Survey (VFTS) 682. This VMS, visible in the Tarantula Nebula (30 
Doradus) of the Large Magellanic Cloud (LMC), has an inferred present-day mass of ~ 150Mq and 
is located at a projected distance o f ~ 30 pc in the North-Ea st from the central massive star cluster 
R136 of the 30 Doradus complex (IBestenlehner et al.ll201ll) . From radial velocity measurements 



and its projected separation from R136, IBestenlehner et al.l (|201ll ) inferred that the true velocity 
of VFTS 682 should be 40 km s"^ if it is indeed a runaway from R136. On the other hand, the 
apparent isolation of VFTS 682 f rom any star cluster, its moderate peculiar motion and the non- 
detection of a bow shock (but see iGvaramadze et al.ll2010l) might also indicate that this VMS has 



formed isolated and is unrelated to R136, as IBestenlehner et al.l (j201ll ) argue. The origin of VFTS 
682 is therefore currently unclear, i.e., whether it is a "slow runaway" from R136 or has formed 
alone, outside any clusters. Thi s, in turn, freshens up the open question whether massive stars 



form only in dense environments (jBonnell et al.ll2004l) or whether isola ted massive star formation is 
possible through pure circumstellar disk-accretion (iKuiper et al.ll2010l ) . Notably, the stellar content 
of R136 itself presents a challenge to the current understanding of star-formation phenomenology 
as it contains several single-star members with inferred initial masses upto ^ 32OM0, i.e., well 
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above the wi dely accepted 1 50Mp canonical upper- limit of the stellar IMF (jCrowther et al.l 12010 ) 



As coined bv lKroupa et al.l (|201ll ). the stellar population of R136 is therefore "super-saturated" by 



containing super-canonical stars. 

Because of the important implications in massive star formation theory, it is essential to deter- 
mine the likeliness of a VFTS 682-like runaway from R136. To that end, we compute the evolution 
of model star clusters with properties resembling R136 using the direct N-body integration method 
and focus on the kinematics of the massive stars ejected from them. We find that the ejection of 
VMSs, within a few Myr, with kinematic properties similar to VFTS 682 is quite plausible: the 
latter star is then quite possibly a runaway from R136 and the substantially controversial assump- 
tion of isolated massive star formation is unnecessary to explain its existence. These computations 
also demonstrate that super-canonical stars appear naturally from merged binaries. 

In Sec. [21 we describe our model clusters (Sec. 12. ip and the method of our computations 
(Sec. 12. 2p . In Sec. [3] we present our results focussing on VMS runaway members that resemble 
VFTS 682 (Sec. 13. ip and also on 30 Dor 016-like runaways (Sec. 13. 2p . We also discuss how the 
ejected fraction of stars depends on the stellar mass (Sec. 13. 3p . We conclude this paper in Sec. [H 
with a discussion and highlighting possible future improvements. 



2. Computations 

2.1. Initial conditions 

As our primary objective is to study the runaway massive stars from a R136-l ike cluster, we 
begin our N-body computations with star clusters modelled as Plummer spheres (jKroupal 120081 ) 
having properties conforming with the observed parameters of R136. The present-day parameters 
of R136 still remain ambiguous. We take the total initial mass of each of our model Plummer 
spher es to be Mci{0) ~ IO^Mq which is an upper limit of the mass of R136 ( Crowther et al. 
20ld ). As for its size, a va riety of half-mass radii /core radii has so far b e en assigned to this 



clust e r by different author s (IHunter et al 



1995 



20031 : lAndersen et al.l l2009l : iPortegies Zwart et al.l |2010|) . We take the initial half- mass radii of 



Portegies Zwart et al.l l2002l : iMackev Sz Gilmore 



our models to be r/j(0) 0.8 pc. The core radii of our models t hen turn out to be < 0.3 pc 
throughout their evolutions, 0.3 pc being an observed upper limit (iMackey fc Gilmord l2003l ) . The 
core density is > 1.3 x 10^ stars pc~^. 



^The initial masses of the clusters' member stars are drawn from the canonical IMF (jKroupa 

200ll ) in the range O.O8M0 < < ISOMq. This IMF is a two-part power-law with indices ai = 1.3 
for 0.08Mq < rus < 0.5Mq and 02 = 2.3 for more massive stars. Throughout this paper, we denote 
the (instantaneous) mass of a stellar member by in general and M is exclusively reserved to 
denote the mass of a runaway (see Sec. [3]) member. The metallicity of the stars are taken to be 
low: Z = 0.5Zq which is appropriate for R136. 
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Most stars form in binaries or higher-order multiple systems. Furthermore, c lose binary-binary 
encounters play an important role in the dynamical ejection of stars from a cluster ([Leonard &: Duncan 
19881 . 119901 ) ■ However, at the current state-of-the-art, binaries are computationally the most time- 
consuming members in a direct N-body calculation of a star cluster because of the particular 
numerical algorithms involved (see Sec. 12. 2p . Thus, in order to achieve a reasonable calculation- 
speed, we initially arrange all the stars more massive than 5Mq in binaries and all the less massive 
members are kept single. As our primary aim is to probe the massive runaway stars, which can be 
efficiently ejected only via close encounters with massive, hard binaries (see Sec. 13. only massive 
binaries are important in this study. 

To create the initial binary population, we first sort the > 5Mq stars with decreasing mass 
and then pair them in order, so that the resulting binaries have their mass-ratios close to unity. It 
is already understood that a mass-ratio distri bution biased towards unity is required in order to 
get a significant ejection rate of massive stars (IClarke fc Pringldll992ll. Furt hermore, observations 



indicate that 0-stars favor OB stars as their companions (jSana &: Evansll201ll ). The orbital periods, 
P, of the binaries with primary masses > 2OM0 are chosen from a uniform distribution between 
0.5 < log^o P < 4 (Opik's law; P in days) since the 0-type sta rs are preferent i ally f ound in short- 
period systems. This range of P is similar to that obtained bv lSana &: EvansI (j201ll ). 



For 5Mq < lUg < 2OM0, equation (8) of iKroupal ()1995bl ) is chosen as the period distribution, 
which is given by 



P, birth 



V 



(logio P 



logio -Pmin) 



6 + (logio P - logio Prr 



\2 ' 



with 7? = 2.5, S = 45 and which cove rs a much wider range of period between 1.0 < logio P < ^-^^ 



(jKroupalll995a|jbl: iMarks et al.ll201ll ). Here, /p,birth is the "birth period distribution" of binaries as 
obtained bv lKroupal (|l995al ) through "inverse dynamical population synthesis". It is the primordial 
binary period distrib ution unmodifi ed by the dynamical destructions of these binaries (i.e., without 
a depletion function; lKroupalll995al ) or by the mutual interactio ns between the binary members in 



their pre-main-sequence stage (i.e., without an "eigenevolution" ; lKroupalll995bl ). The modification 



of a primordial binary's orbital period, mass-ratio and eccentricity following its eigenevolution is 
known only for low mass stars (jKroupa &: Petr-Gotzena I2OIII : iMarks et al.l l201lh . We choose the 



orbital eccentricities e following a thermal distribution, i.e., /(e) = 2e (jSpitzerll 1 987l : lKroupall2008l ) . 
Fig. [1] shows the resulting distribution of binary binding energies as a function of the primary mass. 
The two distinct energy distributions across rus = 20Mq are clearly visible; the binaries with 
nis > 20Mq primaries are harder than those with < 2OM0 primaries. 

The initial systems are completely mass-segreg ated (|SDitzeJl987l ) using a method bv lBaumgardt et al, 



(|2008l ) so that the heaviest stars are initially concentrated within the clusters' cores. Such an initial 
condition mir nics primordial mass segregation that is inferred to be true for several Galactic glob- 
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ular clusters (jBaumgardt et al 



) and young star clusters (jLittlefair et al 



2008; 



Marks fc Kroupa 



2003 



2010; 



Zonoozi et al. 



Chen et al 



2011 



Hasan &: Hasan 



20071 ). Notably, the initial setup 
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employed here (all massive stars in binaries, fully mass-segregated cluster) are the very first ones 
ever computed under such extreme conditions. 



2.2. N-body integration, stellar evolution and stellar hydrodynamics 



The initial models, constructed as above, are dy namically evolved using the state-of-the-art 
direct N-body integrator "NB0DY6" (jAarsethI 120031 ) which takes advantage of the remarkable 
hardware accelerated computing capacity of a Graphical Processing Unit (GPU) while integrating 
the centers of masses. Al gorithmically, the most important feature of NB0DY6 is that it applies 
regularization techniques (jAarsethI |2003| ) to resolve close encounters, making them highly accurate 
and therefore uniquely suitable for this work. However, when there are a significant number of 
primordial binaries, their regualrized orbits can currently be integrated only on the much slower 
host workstation processors, which bottlenecks the GPU's hardware acceleration significantly. No 
external tidal field is applied during the N-body integration as the gravitational field structure of 
LMC is currently unclear. 

In addition to integrating t he point-mass mot ion, NB0DY6 also employs analytical but well- 
tested stellar evolution recipes (jHurlev et al.ll2000l ) to evolv e the individual s tars. These prescrip- 
tions are based on model stellar evolution tracks computed bv lPols et al.l (119981). The wind mass loss 
of all massive stars are taken into account using the empirical formula of lNieuwenhuizen de Jager 
(jl990l ). while they are on the main sequence. The mass loss rates are, of course, appro priately modi- 



fied o n the giant branches and other evolved phases for stars of all masses (se e Sec. 7 of iHurlev et al 



2000). The code incorporates the physics of stellar binary evolution as well (iHurley et al.ll2002l ). It 
also includes detailed models for tidal interactions between stars and recip es for the outcomes of 
mergers between different types of stars and stellar remnants (see Table 2 of iHurley et al.ll2002l for 
a summary). 

Since mergers among main-sequence (hereafter MS) stars have substantial consequences in our 
results (Sec. [3l ), we summarize he re the treatment of a merger betwee n two MS stars in NB0DY6 
(which follows IHurley et al.ll2002l schemes; see also IHurley et al.ll2005l ). When two MS stars collide 
(a collision between two single MS stars or between two MS components in a binary due to eccen- 
tricity induced by close encounters and/or due to encounter hardening), it is assumed that (a) the 
merged product is a MS star with the stellar material completely mixed and (b) no mass is lost 
from the system during the merger process. The no mass-loss assumption is based on the results 
of SPH simulations of MS-MS mergers {e.g., ISills et al.l [2OO1I 1 but such calculations yield only a 
limited amount of mixing. The rejuvenated age of the merged MS star is determined depending 
on the amount of unburnt hydrogen fuel gained by the hydrogen-burning core as a result of the 
mixing. In case of a mass transfer across a MS-MS binary, distinction is made between the cases 
when the original accretor MS star has a radiative or a convective core. For a convective core, the 
core grows with the gain of mass and mixes with the unburnt hydrogen fuel so that the accreting 
MS star appears younger. For the case of a radiative core, the fraction of the hydrogen burnt in 
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the hydrogen-burning core remains nearly unaffected by the gain of mass so that the effective age 
of the MS star decreases. The effective age is determined so as to keep the elapsed fraction of its 
MS lifetime unchanged. 

The initial cluster mass of Mci{0) ~ 10^ Mq in our models corresponds to A'^(O) = 170667 stars. 
To our knowledge, direct N-body computations with such a large number of stars, where the clusters 
are fully mass-segregated and all the massive stars are in binaries, are being reported for the first 
time. We evolve 4 initial models with the above A^(0), generated using diff erent random number 
seeds, until ^ 3 Myr. We ta ke this age as an upper limit of the age of R136 (jCrowther et al.ll20ld : 



Portedes Zwart et al.ll20ld ). We do ah the computations on "NVIDIA 480 GTX" GPU platforms. 



3. Results 

From the above computations, we trace the bodies that are ejected from the clusters during 
their evolution (within ~ 3 Myr). We consider a single-star/binary/multiplet to be a runaway 
member from its host cluster if it is found moving away from the cluster beyond i? > 10 pc 
distance from the cluster's center of density. Although our primary focus is on the runaway VMSs, 
we consider the whole mass spectrum of ejected stars as well. 

Fig. [2] shows the projected snapshots of the runaways with masses M > 3Mq, combined from 
the 4 computations, at t = 1 Myr and 3 Myr evolutionary times. It can be seen that by t = 3 
Myr there are a significant number of fast runaway VMSs with total (3-dimensional) velocities upto 
~ 300 km and a few fast VMSs are already present at t = 1 Myr. All these runaways are on 
their MSs and hence are OB stars. We note that the vast majority of the massive ejected members 
are single stars — only 2 of the massive ejecta from our computations are found in hard binaries. 
We shall discuss the multiplicity properties of the ejected stellar population in detail in a future 
paper. 

It is worthwhile to note that although our adapted canonical IMF has a 15OM0 upper limit, 
our models yield single-star runaways with masses upto ~ 250Mq within t < 3 Myr, considerably 
exceeding the widely accepted 150Mq limit. A few of such members are also found bound to each 
of our model clusters within the same evolutionary perio d. These VMSs form when massive bin ary 



components merge as a result of encounter hardening (jHeggid Il975l : iBaneriee Sz GhoshI 120061 ) of 
these binaries and/or eccentricities induced to them by close encounters. It is the adapted relatively 
small and narrow orbital period distribution of the massive binaries (m^ > 2OM0), conforming with 
the observed properties of 0-type stellar binaries (see Sec. 12. ip . that makes such mergers probable, 
leading to the formation of single stars with > 150Mq. On the other hand, it is these hard, 
massive binaries that can efficiently drive the VMS runaways (see Sec. 13. ip . In other words, the 
possibility of the presence of runaway VMSs naturally leads to the formation of VMSs (bound 
or ejected), in the course of the evolution of the cluster, with masses significantly exceeding the 
canonical upper limit, even though the cluster's IMF is intrinsically canonical. Therefore, the 
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observed super-canonical VMSs in R136, with inferred individual initial masses upto ~ 320Mq, do 
no t necessarily imply an IMF with an upper-limit greatly exceeding the canonical limit as argued 
by ICrowther et alj (|20ld ): the > 150Mq single-star members can as well be accounted for as 



recent massive binary merger products. 

Each of our computed models yields a fe w of such sup er-canonical single-star members. Because 
of their substantial gravitational focussing (|Spitzeiill987l ). they are involved in very frequent close 
encounters with other massive binaries, thereby being vulnerable to being ejected from the cluster. 
Three of our 4 computed clusters have produced super-canonical single-star runaways with M > 
I5OM0 and the most massive merger product has been ejected from 2 of the models. 



3.1. VFTS 682: a very likely runaway from R136 

To have an understanding of the spectrum of velocities with which the stars run away, we plot 
their 3-dimensional velocities V as a, function of their instantaneous masses M bX t = 1 Myr and 
3 Myr as shown in Fig. [3l The data points from all the computations have been compiled in this 
figure, which are distinguished by the different symbols. We note two trends in this plot, viz., (a) 
there is a lower boundary of the scatter in V that increases moderately with M and (b) the upper 
boundary of this scatter is independent of M. These trends are vividly apparent from the plot 
at t = 3 Myr (Fig. [3l lower panel) where there is a significantly larger number of runaways. No 
well-defined boundaries can, of course, be drawn in this plot and the above mentioned boundaries 
are meant to indicate net trends. 

To understand these trends, we first note that an ejected single star of mass M and velocity 
V can appear in two main w ays, viz., (i) it is ejected due to a kinetic energy (K.E.) boost in a 



close, super-elastic encounter (|Heggiel Il975l ) with a hard binary where no exchange of members 
has occurred (a fiyby encounter) and (ii) it is ejected from a hard binary being replaced by a 
more massive star (an exchange encounter). In both cases, the K.E. of ejection is of the order of 
the binary's binding energy, on average. The trend of the lower boundary can be understood by 
noting that a single star is likely to experience a flyby in a strong encounter with a hard binary 
(most of them have nearly equal mass components in our model) only if the binary components 
are more massive than the intruder, since otherwise the latter is likely to get trapped in the binary 
by exchanging with one of its (lighter) members. In other words, to recoil a mass M in a flyby 
encounter, one needs a binary with components of mass rus ^ M . Similarly, a star of mass M can 
preferentially be ejected from a binary only by being replaced by an intruder of mass > M. Since 
the ejected star, in both cases, w ill typically have a K.E. of the order of the binary's binding energy 



( Heggielll975l : iHeggie et al.lll996l ). we have MV"^ ~ ml/a, where a is the binary's semi-major-axis. 
Hence the condition > M implies V"^ > M/a; the smallest possible ejection velocity increases 
with M, for a given a. Hence, the widest (hard) binaries in our model give rise to the lower 
boundary of the scatter in Fig. [3] (lower panel). Indeed, the eye-estimated lower limit of the data 
for the lower mass ejecta can be well represented by a curve V oc M^/^ shown by the pink line 
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segment in Fig. [3] (lower panel). This is also true for the massive end; their lower boundary can 
be represented by another curve (the green line segment) but with the same proportionality. The 
discontinuity is expected due to that in the binary period distribution (c./. Sec. 12. H Fig. [T|). 

On the other hand, the upper boundary of the plot is formed by the ejecta generated by 
binaries having the highest binding energies. Fig. [T] shows that the hardest binaries are also the 
most massive ones, as a result of our chosen primordial binary distribution (see Sec. 12. ip . They 
can therefore eject stars over the entire mass range with K.E. of the order of their binding energy 
forming the upper boundary that is independent of M. Note that as an order of magnitude estimate, 
the highest binding energies are ~ 10^ - W'^Mq pc^ Myr-2 which yield V ^ lO^'^ - 10^ km s'^. 
This indeed agrees with the highest Vs in Fig. [3] in terms of order of magnitude. 

In Fig. [3] (lower panel) it can be seen that the median of the ejection velocities (thick black 
line) tend to F ~ 50 km s"^ for the most massive runaways and the corresponding quartiles (thin 
black lines) span between about 40 — 60 km s~^. Taking into account the outliers with V smaller 
than the first quartile in the massive end of the M — V plot, these ejection velocities are similar to 



the inferred true velocity of VFTS 682 (jBestenlehner et al.ll201ll ). This result, therefore, already 



strongly indicates that VFTS 682 is a typical runaway VMS from R136. The highest V reached in 
our computations, however, exceeds 300 km as can be read from Fig. [3l 

We further inspect that all of our computations eject VMSs with kinematic properties agreeing 
fairly with those of VFTS 682, viz., F 40 km s^\ i? 30 pc and M ^ ISOMq. Such instances 
are shown in Table [TJ These results dictate that the VFTS 682 is an expected, very likely runaway 
VMS from R136. 

As additional information to the reader and for convenience of comparison with observations, 
we provide the data corresponding to Fig. [3] in the form of an online table (Table [2]) • In addition to 
the 3-dimensional velocities and the (instantaneous) masses of all the runaways from all of our com- 
putations, we also provide the corresponding line-of-sight and tangential velocities and projected 
positions. For the aid of the reader in doing further analyses, we also provide the corresponding 
components of the positions and the velocities in an additional online table (Table [3]). 



3.2. Another runaway: 30 Dor 016 



Another notable VMS in the 30 Doradus is 30 Dor 016 which is also thought to be a runaway 
from the central R136 cluster (jEvans et al.ll20ld ). This star, located at ~ 120 pc projected distance 
from R136, has a radial velocity of ~ 85 km and its mass is inferred to be ~ 90Mq (see 
Evans et al.ll20ld and references therein). The possibility of #016 having a close companion being 
quite unlikely as these authors conclude based on VFTS multi-ep och spectroscopy, it is very likely 
a runaway from R136. Given that this VMS is about 1 Myr old (lEvans et al.ll20ld ). its projected 
distance implies a minimum of ~ 120 km s^^ tangential velocity. Combining this with its radial 
velocity, #016 has at least ~ 150 km s^^ 3-dimensional velocity. Two of our computations are 
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indeed found to eject runaways with similar properties which are shown in Table HI 



3.3. Mass-dependence of runaway stars 



It is very interesting to study the mass (or spectral class) dependence of the runaway members 
as it directly points to the pure dynamical nature of the ejection process. Although the number 
of stars bound to the cluster decreases strongly with increasing stellar mass, more massive stars 
and binaries are more centrally concentrated due to mass segregation and hence interact with each 
other more efficiently and frequently. This causes the runaway fraction of stars (defined as the 
ratio of the number of stars, within a bin around mass M, moving away from the cluster at i? > 10 
pc to the total number of stars in the whole system, i.e., including all the bound and the ejected 
members, within the same mass bin) to increase with the stellar mass as shown in Fig. 01 where the 
outcomes of all the computations are overplotted at t = 1 Myr and 3 Myr. In Fig. [H the ejected 
fraction g{M) increases considerably with M for M > bM^ but remains nearly flat and small for 
lower M. This is probably an artifact of our initial binary population where only the stars with 
rus > 5Mq are in binaries and therefore are efficient in ejecting only those members which are more 
massive than 5Mq. Current technology does not allow us to perform direct N-body integrations of 
R136-type clusters in which all stars are initially in binaries. S uch computations are available only 
for moderate mass clusters (jKroupalll998l : iKroupa et al.ll200lh . 



4. Discussions and outlook 

In the present work, we have computed the dynamical evolution of model clusters whose 
structure conform with the observed global properties of R136 — the central massive cluster in the 
30 Dor complex of the LMC, using the direct N-body integration method. The evolution of the 
i ndividual stars, chosen i nitially from the canonical IMF with the standard 150Mq upper cutoff 



(jWeidner Sz Kroupal l2004l ) . has also been incorporated and as well the evolution of the individual 
binaries. We focus on the ejection of massive stars from our model clusters which is a process that 
depends crucially on the properties of the primordial binaries in the cluster. For computational 
simplicity, we have taken only the stars with initial masses > 5Mq to be in binaries (see Sec. l2.ip 
which are the only ones that are efficient in ejecting the massive stars. Hence, the properties of 
the massive runaways are not expected to be affected significantly by the absence of lower mass 
binaries. 

Recent observations of the R136 and the 30 Dor region have raised fundamental questions 
regarding massive star formation mechanisms. In particular, the apparently isolated and relatively 
slow-moving single VM S VFTS 682 has raised s uspicion that it might be an instance of isolated 



massive star formation ( iBestenlehner et al.ll201ll : see Sec. [T|). Additionally, the presence of single 
VMSs in R136 with inferred initial masses upto « 32OM0 has questioned the canonical 150Mq 
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upper limit of the IMF (jCrowther et alJboid ) 



The most important outcome of our computations is the confirmation that a "slow runaway" , 
with a 3-dimensional velocity similar to that of VFTS 682, is in fact the most probable type of 
ejected VMS from a R136-like cluster (c./. Fig. [3l Sec. I3.ip . In fact, all of our computed models 
yield one or more runaways with kinematic properties agreeing fairly with those of VFTS 682 (c./. 
Tabled]). Given such a likeliness of a VFTS 682- type runaway from R136, this apparently isolated 
star clearly does not imply isolated massive star formation and it is very likely a former member 
of R136. 

Furthermore, as explained in Sec. [3l massive, close binaries are necessary to dynamically eject 
VMSs from star clusters, which, in turn, are susceptible to merge due to their hardening and/or 
eccentricity-pumping, by the frequent close encounters that they receive. As our computations 
show (see Sec. [3]), such massive binary mergers can easily produce single stars, within a few Myr, 
with masses well exceeding the 150Mq upper limit, even if the cluster begins with the canonical 
upper limit. Our 4 computations have produced upto ~ 250Mq members and merger products 
upto ~ 300Mq are possible if the most massive binaries merge. Therefore, it doesn't seem to 
be a surprise that R136 has upto ~ 32OM0 single-star mer nbers, given the larg e uncertainties 
in the stellar evolution models used to infer the masses (see ICrowther et al.l l201Cl and references 
therein) and the presence of these super-canonical VMSs (or the presence of a super-saturated 
mass function) is not an indication that R136 was born with an IMF having an upper limit that 
substantially exceeds the widely accepted ISOMq limit. These VMSs can as well be massive binary 
mergers. 

It is, of course, important to note that there are aspects in our assumptions and initial condi- 
tions that favor ejection of VMSs and formation of massive merger products. First, the assumption 
of no mass-loss during MS-MS mergers (see Sec. 12. 2p is an idealization. Al though SPH computa- 



tions indicate that no mass is lost during a MS-MS collision (jSills et al.ll200ll ) , the studied collisions 



are typically between low mass stars and such mass conservation is not necessarily true for massive 
MS-MS mergers. Furthermore, the merged MS star can spin rapidly, even beyond break-up, in 
which case the net mass gain would be much smaller or nil. Second, the assumption of a com- 
plete mixing (see Sec. 12. 2p is also an idealization; it is possible that the denser helium-rich core(s) 
preferably sink to the center of the merger product resulting only a partial mixing of unburnt hy- 
drogen with the helium core. This would lead to a partial refuelling or rejuvenation of the merged 
MS stars's core and hence a smaller lifetime. In other words, the treatment of complete mixing 
maximizes the lifetime of the merged star and hence its probability of getting dynamically ejected 
and/or being observationally detected. Although idealized, these conditions are those which can 
currently be adapted at best i n a direct N-bo dy calculation; other wid ely-used direct N-body codes 
such as the "NB0DY6++" dSpurze J flQQfll ) and the "STARLAB" jPortegies Zwart et al.l \m)^ ) 
also utilize similar synthetic stellar evolution and merger recipes. 



Third, the VMS ejection and the massive merger formation are also facilitated by our adapted 
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initial or primordial mass segregation (see Sec. I2.ip . While the motivation for choosing this condi- 
tion is indeed to maximize these effects, primordial mass segregation is inferred to be true for se veral 



Galactic globular c lusters (IBaumgardt et al 



2008; 



Marks Krqupa 



2010 



Strader et al.ll201lh and 



open clusters {e.g., Bonnell Davie J 1998 : Littlefair et al. 2003 : Ichen et al. 2007 ). An important 



outlook would be to study the effects of a varying degree of initial mass segregation. 

Another ingredient that causes ejection of runaway VMSs and the presence of nis > 150Mq 
VMSs in our computations likely, is the presence of relatively close primordial binaries beyond 
nis > 20Mq (see Sec. l2.ip . However, such a binary distribution is in accordance to the observational 
fact that 0-stars are generally found in close binaries. Notably, the range of the 0-star binaries' 
initial orbital periods in our models is taken to be similar to the observed one (jSana &: Evans 
). In this context, recall that we adapt a P-distribution that follows a single Opik's law and 



2011 



also take the binary mass-ratios (q) close to unity (see Sec. 12. ip . However, the P-distribution of 
0-stars in 0-rich clusters is observed to follow a bi-uniform distribution in log^^g ^ with the break 
at P ~ 10 days; there is a significant overabundanc e of binaries (compri sing 50%-60% of the 0-star 



binary population) for P < 10 days (see Eqn. 5.2 of lSana &: Evansll201in. Furthermore, th eir mass- 



ratios are found to be uniformly distributed between 0.2 < q < 1.0 (jSana &: Evansll201lh . A very 



important development would be to incorpo rate such more r ealisti c distributions of 0-star binaries' 
orbital periods and mass-ratios. In fact, the lSana &: EvansI (|201ll ) P-distribution would even more 
strongly favor the ejection of VMSs and the formation of super-saturated stars than the presently 
adapted distribution due to the overabundance of P < 10 days binaries in the former distribution. 
This would, of course, be at least partially counter-balanced by the g-distribution; the smallest q 
gives a binary with a binding energy smaller by a factor of a few for the same primary-mass and 
semi- major- axis. It is, therefore, an essential outlook to study the net outcome of such effects. 
However, note that the gross conclusions, as obtained from the present models, are unlikely to get 
affected by these details since the range of P and hence the range of binary binding energies are 
indeed taken to be similar. 

Finally, note that the c hosen masses of the in itial Plummer clusters are the upper mass limit 
of that of R136 (?^ lO^M rr^; ICrowther etaP boiol ). A lower mass limit of R136 is pa 5 x IO^Mq 



(jWeidner &: Kroupall2004l ). This being only by a factor of two smaller than the upper limit, we do 
not expect that the results presented here would be affected significantly had we chosen the lower 
mass limit as the total initial cluster mass. 

A corollary of our above two results is that runaway VMSs and super-canonical VMSs (with 
masses exceeding the canonical upper-limit) would coexist, in general. Of course, as explained in 
Sec. El the latter type of VMSs can easily be ejected from their host clusters and may be found only 
as runaways. It would be worthwhile to search for other massive (and young enough) star clusters, 
in addition to R136, where such super-canonical VMSs exist and the clusters can be related to VMS 
runaways at the same time. Incidentally, the VMS ( Pup (~ 70Mq), whose kinematics indicate 
that it is a runaway from the Vela R2 associatio n, has recently been inferred to be a merger product 
of at least two massive stars ( VanbeverenI 120111 ) . 
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The main drawback of this study is the incompleteness of the primordial binary population. As 
already stressed earlier, the absence of initial binaries for < SM© is anyway unlikely to signifi- 
cantly influence the ejected stellar population for M > 5Mq which justifies such a truncation for our 
purposes. The truncated primordial binary population, of course, has effects on less massive ejecta 
which are pointed out in Sec. [3l The purpose of eliminating the lower mass binaries is, of course, 
computational ease. The presence of binaries substantially diminishes the speed of direct N-body 
computations (even on GPU platforms; see Sec. 12. 2p and such computations with a full spectrum of 
primordial binaries (i.e., 100% primordial binary fraction) of models of the mass that we use in this 
study (i.e.. A'' ^ 170000; see Sec. l2.2|) are currently prohibitive even with the fastest available GPUs 
and workstation processors. Nevertheless, we wish to emphasize again that our computations in this 
study are the largest sized direct N-body computations of initially fully mass-segregated clusters in 
which all massive stars are in binaries rep orted so far, thanks to the remark able boost of computing 



speed provided by the recent G PUs (see INVIDIA GEFORCE Homepage! ) and to the efficient nu- 



merical algorithms of NB0DY6 (jAarsethll2003l ). These calculations represent the first-time attempt 



to directly evolve a real-sized massive star cluster. 

An immediate improvement over our work is, of course, to explore a full primordial binary 
spectrum and the effect of their varied distribution functions which, we expect, will become plausible 
in the near future given the continuing remarkable algorithmic improvements of NB0DY6 and the 
availability of increasingly faster GPUs. A more accurate structural modelling of R136 is also 
pending once its dimensions are more clearly settled (see Sec. 12. ip . 

In conclusion, our calculations imply that the observed VMSs outside R136 and also its bound 
VMS members are fully consistent with the standard paradigm of massive star formation exclusively 
in dense environments and following the canonical IMF. Hence, the observations discussed here 
cannot be considered as instances violating these regimes. 

We are thankful to Sverre Aarseth of the Institute of Astronomy, Cambridge, U.K., for his 
effort in the continuing remarkable improvements of NB0DY6 and his gracious help during our 
computations. We are also thankful to Keigo Nitadori of RIKEN, Japan, for his developments that 
resulted in the exemplary hardware-acceleration of NB0DY6. We thank the anonymous referee 
whose suggestions have improved the quality of the manuscript substantially. 
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Fig. 1. — The binding energies of the initial binaries in our computed models as a function of their 
primary masses, where the two groups of binaries are clearly distinguishable across 20Mq. The 
binaries from the 4 initial models (see Sec. 12. 2p are superimposed which are distinguished by the 
different symbols. Binaries with primary masses > 20Mq are assigned significantly smaller and 
narrowly distributed orbital periods than the other group (see Sec. 12. ip making the former ones 
significa ntly harder. The ran ge of binding energies compares well with that of observed 0-star 
binaries (Sana fc Evans 2011). 
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Fig. 2. — Snapshots of runaway stars (a single filled circle) with M > SMq at t = 1 Myr (top) 
and 3 Myr (bottom) evolutionary times as obtained from our computations. The snapshots from 
the 4 computations are superimposed in each panel. The stars are colour-coded according to their 
masses at age the values in the colour code (colour-bar on the right hand side) being in Mq. The 
direction of the arrow that originates from each star gives its direction of motion in the plane of the 
snapshot while its length is proportional to the 3-dimensional velocity of the star, the scale being 
shown at the upper-left corner of each panel. It can be seen that there are a significant number of 
fast runaway VMSs by t = 3 Myr with 3-dimensional velocities reaching upto ~ 300 km s~^ and a 
few fast VMSs are already present at t = 1 Myr (also, see Fig. [3]). See text for details. 
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Fig. 3. — The 3-dimensional velocity V vs. the mass M of all the runaway single stars (filled 
symbols) at t = 1 Myr (top) and 3 Myr (bottom) as obtained from our calculations. Different 
symbols represent the outcomes from different computations. The thick black line in the bottom 
panel is the median of the scatter along the V axis and the thin black lines are the corresponding 
first and the third quartiles. In computing these percentiles, the data-points are divided into 15 
equal bins in the logarithmic scale over the range O.OSM© < M < 3OOM0. At the massive end, the 
scatter lies within the quartile range of about 40 - 60 km (bottom panel) in striking similarity 
with the inferred true velocity of VFTS 682. The lower boundary to the majority of the data-points 
can be well represented by the power-law V oc M^/^, which are the slim pink and the green line 
segments (bottom panel). The percentiles and the lower boundaries are constructed only for the 
bottom panel (t = 3 Myr) as the data in the top panel {t = 1 Myr) is much sparser. See text for 
further details. 
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Fig. 4. — The runaway or ejected fraction of stars g{M) as a function of stellar mass M at t = 1 
Myr (top) and 3 Myr (bottom). The masses are divided into 25 equal bins in the logarithmic 
scale over the range O.OSM© < M < SOOMq. The outcomes from all the computations, which 
are distinguished by the different symbols, are overplotted in each panel. For M > SMq, g{M) 
increases considerably and becomes nearly unity for the most massive bin. This indicates that the 
dynamical ejection process becomes increasingly important with increasing stellar mass such that 
the most massive members are nearly all ejected. The flatness of g{M) for M < 5Mq may be an 
artifact of the primordial binary population adapted in the computations. See text for details. 
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Table 1. A list of runaway single stars with properties fairly close to those of VFTS 682, as 
obtained from our computations. The columns are as follows: Col. (1): model number, Col. (2): 

evolutionary time t at which the runaway is detected. Col. (3): mass M of the runaway at t, 
Col. (4): its distance R from the cluster's center of density and Col. (5): its 3-dimensional 
velocity V. M, R and V of these runaway stars agree fairly with those e stimated observationally 



for VFTS 682 (last line: Bestenlehner et al.ll2011 



model number 


time t (Myr) 


mass M (Mq) 


distance R (pc) 


velocity V (km s ^) 


1 


2.8 


256.4 


31.9 


27.5 




3.2 


135.9 


26.6 


34.8 


2 


2.6 


126.4 


27.7 


45.7 




2.6 


125.9 


29.9 


49.4 


3 


2.6 


106.9 


45.7 


27.3 


4 


1.9 


169.1 


29.3 


29.0 




1.9 


116.9 


35.2 


32.8 


VFTS 682 


< 3.0 


150.0 


^ 30.0 


^ 40.0 
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Table 2. Data for all the runaway stars as obtained from all of our computations at 1 Myr and 3 
Myr. The columns are as follows: Col. (1): mass M of the runaway at 1 Myr/3 Myr, Col. (2): its 
tangential velocity vt (all velocities are relative to the cluster's center of mass)^ , Col. (3): 
line-of-sight or radial velocity frad*^ i Col. (4): 3-dimensional velocity V , Col. (5): its projected 
distance r from the cluster's center of density, Col. (6): 1 ^ single, 2 binary and Col. (7): 

identity of the star /. 



M (Mq) 


Vt (km s "'") 


Ui-ad (km s 1) 


V (km s-i) 


r (pc) 


single/binary 


/ 


5.9 


60.0 


-25.9 


65.3 


15.3 


1 


1626 


8.6 


10.4 


10.7 


15.0 


9.1 


1 


960 


0.2 


8.4 


8.8 


12.2 


7.9 


1 


4491 


0.4 


15.8 


44.8 


47.5 


5.8 


1 


6332 


0.3 


6.8 


3.3 


7.6 


9.6 


1 


8042 


40.5 


6.5 


-8.2 


10.4 


19.5 


2 


170968 


94.0 


50.3 


209.9 


215.8 


83.7 


1 


10 


236.3 


31.3 


-72.8 


79.2 


9.4 


1 


5 



^The reference frame is chosen fixed and originated at the cluster's center of mass. Its axes are taken to 
be oriented arbitrarily as there is no preferred directionality due to the absence of an external field. All the 
projections are taken on the x — y plane. 

'^The velocity component normal to the plane of projection or the z velocity component is simply taken 
to be the radial velocity due to the arbitrariness in the choice of the coordinate axes. The corrections due 
to the spatial extent is negligible for the distance of R136. 

Note. — The full Table [5] can be accessed online in the electronic edition of the journal (see 
http:/ /tinyurl.com/3kdy64g| for a text version). Examples of its data are shown here to indicate its for- 
mat. 
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Table 3. Data for all the runaway stars as obtained from all of our computations at 1 Myr and 3 
Myr. The columns are as follows: Col. (1): mass M of the runaway at 1 Myr/3 Myr, Col. (2): its 

identity /, Col. (3-5): components of its position relative to the cluster's center of density'^ , 
Col. (6-8): components of its velocity relative to the cluster's center of mass, Col. (9) 1 =^ single, 

2 =^ binary. 



M (Mo) 


/ 


X (pc) 


y (pc) 


z (pc) 


Vx (km s ^) 


Vy (km s -"-) 


tiz (km s ■"-) 


single/binary 


5.9 


1626 


-12.6 


-8.7 


-6.5 


-49.7 


-33.6 


-25.9 


1 


8.6 


960 


-6.9 


5.9 


9.0 


-8.0 


6.7 


10.7 


1 


0.2 


4491 


4.4 


6.6 


8.8 


3.4 


7.6 


8.8 


1 


0.4 


6332 


-2.2 


5.3 


16.4 


-6.5 


14.5 


44.8 


1 


0.3 


8042 


-6.8 


-6.7 


4.5 


-4.7 


-4.9 


3.3 


1 


40.5 


170968 


-16.1 


-11.1 


-25.2 


-5.4 


-3.5 


-8.2 


2 


107.2 


170726 


10.7 


-55.2 


47.6 


6.3 


-32.7 


28.2 


2 


236.3 


5 


-8.8 


-3.4 


-94.9 


-7.8 


-30.3 


-72.8 


1 



'^The reference frame is chosen fixed and originated at tlie cluster's center of mass. Its axes are taken to be oriented arbitrarily 
as there is no preferred directionality due to the absence of an external field. 

Note. — The full Table |3] can be accessed online in the electronic edition of the journal (see the URL 
|http://tinyurl.com/6fepz6b for a text version). Examples of its data are shown here to indicate its format. 



Ta ble 4. A list of ru naway single stars with M, R and V similar to those of 30 Dor 016 (last 
line: lEvans et al.ll2010l and references therein), as obtained from our computations. The meanings 

of the columns are the same as in Table [TJ 



model number 


time t (Myr) 


mass M (M©) 


distance R (pc) 


velocity V (km s 


3 


3.1 


90.9 


103.3 


144.9 


4 


1.3 


138.6 


122.0 


131.8 


30 Dor 016 


« 1.0 


^ 90 


« 120 


K 150 



